Method for modelling a surface and device for implementing same

ABSTRACT

A method for obtaining a model of a surface including the steps of obtaining measurements of geometrical data concerning specific points on the surface, making a grid of the surface, with the grid passing through said points, memorizing, at an address which is specific to each node of the grid, the coordinates of the node, the number of satellites of the node, information for access to the addresses of said satellites and thereafter to information which relates to them, and geometrical data which may be associated with said node; for each node, defining a local roughness index obtained from a weighted sum of the current coordinates of the node and its satellites, defining the sum of an overall roughness index representing the sum of all the local roughness indices, and of an overall index of the infringement of said geometrical data, iteratively adjusting the coordinates of indefinite nodes, by using at each adjustment the sum of a weighted combination of current node neighbor coordinates and of a combination of geometrical data associated with the node, in order to minimize said sum, and creating a model of the surface on the basis of the adjusted coordinates.

SUMMARY OF THE INVENTION

The present invention is generally concerned with investigations withina three-dimensional body and is particularly concerned with a newprocess of obtaining a representation of a surface within a body on thebasis of a limited set of known geometrical data relating to saidsurface.

Investigations within three-dimensional bodies are of major concern ingeology, geophysics and biology.

In geophysics, for example, it is necessary to obtain as accurate aspossible representations of surfaces situated, for example, at theinterface between two areas of different kinds or with differentproperties, on the basis of data obtained from prospecting, or duringexploitation of an underground resource.

How effectively the three-dimensional body is investigated, especiallyin prospecting for oil or other underground resources, depends on theaccuracy with which this type of surface can be reconstituted andrepresented.

In medicine and biology there are various known processes for obtainingrepresentations of cross-sections through a living body or the like andthe aim is also to reconstitute and represent in three dimensions asurface situated at the interface between two media, for example thecontours of a bodily organ.

There already exist a number of modeling techniques for obtaining arepresentation of a surface within a three-dimensional body that is tobe exploited or treated. Particularly worthy of mention are the Bezierinterpolation method and the spline functions method (for more detailsreference should be had to Geometric Modeling, M. E. Mortenson, Ed. JohnWiley, 1985).

However, all these known techniques are ill suited to handlingheterogeneous data, by which is meant data relating to the coordinatesof points on the surface, to fixed orientations of the surface, togeometrical links between two separate points, etc. Also, thesetechniques are ill-suited to a wide spread of surfaces or to surfaceanomalies such as folds and breaks, and to discontinuous surfaces. To bemore precise, in some cases the functions used in these techniques arenot convergent or have no solution or have no single solution.

Finally, the known techniques are not able to take into account theconcept of the degree of certainty with which some surface data isknown.

The present invention is directed to alleviating the drawbacks of theprior art processes and to proposing a process that can allow fornon-homogeneous and/or highly diverse geometrical data and for theanomalies mentioned above and other anomalies, providing in each case asingle and unambiguous model of the surface.

Another object of the invention is to propose a process that can usegeometrical data and at the same time data as to the degree of certaintyor accuracy of said data.

To this end, the present invention firstly consists in a process formodeling a surface representing for example the interface between twoareas of different kinds or with different properties in athree-dimensional body such as a geological formation or a living body,characterized in that it comprises the steps of:

obtaining by means of measuring apparatus a set of geometrical datarelating to the surface and associated with respective points on saidsurface;

meshing the surface so that all said points are a subset of the nodes ofthe mesh;

storing at a specific memory address for each node of the mesh, thefollowing data:

the coordinates of the node in question,

the number of satellite nodes of the node in question,

data providing access to the specific addresses of said satellite nodesand consequently to the data relating thereto,

if necessary, geometrical data associated with said node in question,

for each node of the mesh, defining a local roughness index derived froma weighted sum of the actual coordinates of the node and of itssatellites,

defining the sum of a global roughness index obtained by summing thelocal roughness indices associated with each node and a global index ofviolation of said geometrical data,

fitting the coordinates of each node for which the precise coordinatesare not known by an iterative method in which for each step of theiteration there are added a weighted combination of the actualcoordinates of the satellite and of the satellites of the satellites ofsaid node and a combination of the geometrical data associated with saidnode, in such a way as to minimize said sum, and

creating a representation of the surface from the fitted coordinates ofeach node.

The invention also concerns a device for modeling a surface representingfor example the interface between two areas of different kinds or withdifferent properties in a three-dimensional body such as a geologicalformation or a living body, characterized in that it comprises:

means for obtaining by means of measuring apparatus a set of geometricaldata relating to the surface and associated with respective points onsaid surface;

means for meshing the surface so that all said points are a subset ofthe nodes of the mesh;

means for storing at a specific memory address for each node of themesh, the following data:

the coordinates of the node in question,

the number of satellite nodes of the node in question,

data providing access to the specific addresses of said satellite nodesand consequently to the data relating thereto,

if necessary, geometrical data associated with said node in question,

calculation means for fitting the coordinates of each node for which theprecise coordinates are not known by an iterative method in which foreach step of the iteration there are added a weighted combination of theactual coordinates of the satellite and of the satellites of thesatellites of said node and a combination of the geometrical dataassociated with said node, in such a way as to minimize a sum of aglobal roughness index obtained by summing local roughness indicesassociated with each node and each derived from a weighted sum of theactual coordinates of the node and of its satellites and a global indexof violation of said geometrical data,

means for creating a representation of the surface from the fittedcoordinates of each node.

These and other objects of the present invention will become morereadily apparent from the detailed description given hereinafter.However, it should be understood that the detailed description andspecific examples, while indicating preferred embodiments of theinvention, are given by way of illustration only, since various changesand modifications within the spirit and scope of the invention willbecome apparent to those skilled in the art from this detaileddescription.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention will emerge more clearly from the followingdetailed description of one preferred embodiment of the invention givenby way of non-limiting example and with reference to the appendeddrawings, in which:

FIG. 1 shows a typical graph incorporating meshing used to model asurface, this graph being defined by the set of all the sides of facetswhose vertices constitute the set Ω described later, and also showsarrows symbolically representing displacements of nodes of the meshapplied interactively by the user to shape the surface;

FIGS. 2(a) through 2(d) show four instances of controlling the shape ofthe surface using four different types of vector constraint; the vectorsΔ.sub.λμ are assumed to be known directly in FIGS. 2(a) through 2(c) andknown because they are orthogonal to a given vector V in FIG. 2(d);

FIG. 3 shows a typical complex geological surface (salt dome) modeledusing the process of the present invention;

FIG. 4 shows a typical complex biological surface (embryo brain) modeledusing the process of the present invention from discrete data obtainedfrom a succession of cross-sections;

FIG. 5 shows the organization of the memory associated with an atom A(k)of the kernel φ_(k) ;

FIG. 6 shows a first example of memorizing an orbit of a given atom inan array with n_(k) items;

FIG. 7 shows a second example of memorizing an orbit of a given atom inn_(k) arrays each of two items;

FIG. 8 shows the memorization of a geometrical datum or constraint;

FIG. 9 shows the method of memorizing a set of geometrical data orconstraints relating to a given atom A(k);

FIG. 10 shows the memorization in an array of all of the atoms relatingto a surface S; and

FIG. 11 shows the chaining of access to said set of atoms;

FIG. 12 shows an apparatus embodiment of the present invention.

Identical or similar elements or parts are designated by the samereference numbers in all the figures.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENT

Reference will be made later to a number of publications whoserespective contents are deemed to be included by way of reference intothis description.

(1) A method of Bivariate interpolation and smooth surface fitting basedon local procedures, H. AKIMA, J. ACM 17.1, 1974;

(2) Shape reconstruction from planar cross-sections, J. D. BOISSONNAT,ACM Transactions on Graphics, vol. No. 2, 1986;

(3) Machine contouring using minimum curvature, I. C. BRIGGS, Geophysics17.1, 1974;

(4) Triangular Bernstein-Bezier patches, G. FARIN, Computer-aidedgeometric design, 1986;

(5) Primitives for the manipulation of general subdivisions and thecomputation of Voronoi diagrams, L. GUIBAS and J. STOLFI, ACMtransactions on Graphis, vol. No. 2, 1985;

(6) Three-dimensional graphic display of disconnected bodies, J. L.MALLET, Mathematical Geology, vol. 20, no. 8, 1988;

(7) Geometrical modeling and geostatistics, J. L. MALLET, Proceedings ofthe Third International Geostatistics Congress, Kluwer AcademicPublishers, 1989;

(8) Discrete Smooth Interpolation, J. L. MALLET, ACM Transactions onGraphics, April issue, 1989;

(9) Geometric Modeling, M. E. MORTENSON, John Wiley, 1985.

Also, the following description of one preferred embodiment of theinvention is followed by section 2, infra explaining a method ofreconstituting a surface derived from the method explained in referencedocument (8) and on which the process of the present invention is based.

1. DESCRIPTION OF A PROCESS IN ACCORDANCE WITH THE INVENTION FORENCODING AND FITTING COMPLEX SURFACES 1.1 Introduction

1.1.1 Nature of the process

An industrial process described hereafter is used to encode and fitcomplex surfaces encountered, for example, when modeling:

interfaces between different domains,

interfaces between geological layers,

interfaces between biological organs,

etc.

The result of this encoding phase may be converted into pictures and/orsolid models (molded from plastics material, for example) by a graphiccomputer or any other device. In geology, for instance, these picturesand/or solid models are used to optimize the prospecting, exploitationand management of earth resources such as:

oil deposits,

mineral deposits,

water deposits,

etc.

The process described hereafter consists of two complementarysub-processes:

1. a sub-process to encode a surface,

2. a sub-process to fit a surface to the data currently available. Ingeology, for example, these data may be:

the exact or approximate location of points in the intersection of agiven surface and a plane,

the exact or approximate location of the intersection of a well and ageological surface corresponding to the interface between two geologicallayers,

the exact or approximate location of the plane tangent to a geologicalsurface corresponding to the interface between two geological layers,

the exact or approximate throw vector of a fault breaking the geologicalsurface corresponding to the interface between two geological layers,

etc.

The encoding process is designed to optimize the fitting process basedon a new version of the DSI method (cf. [8]) especially developed forthe geometrical modeling of surfaces. This new version of the DSImethod, which has not yet been published, is presented in section 2,infra.

1.1.2 Dividing a surface into facets

Let S be the surface to model. As shown in FIG. 1 and as explained in[6] and [7], this surface will be approximated by a set of triangularand/or polygonal facets. These facets are characterized by:

The coordinates of their vertices. These vertices constitute a finiteset of N points {φ₁, . . . ,φN } regularly distributed on S and numberedfrom 1 to N. The coordinates of the k^(th) vertex φ_(k) will be denoted(φ_(k) ^(x),φ_(k) ^(y),φ_(k) ^(z)).

Their edges constituting the links between pairs of vertices(φ_(i),φ_(j)).

The set of edges of the triangles and/or polygons composing the surfaceS constitute a graph G whose nodes are the vertices of the trianglesand/or polygons. For this reason and for the sake of simplicity, "node"is used hereafter as a synonym of "vertex of triangles and/or polygons".

    node↑{vertex of triangles and/or polygons}

1.1.3 Preliminary definitions

Set of satellites Λ(k)

Let φ_(k) be the k^(th) node of a surface S. The "set of satellites" ofφ_(k) is the set Λ(k) Of the nodes φ_(j) different from φ_(k) and suchthat (φ_(j),φ_(k)) iS an edge of triangle or polygon. In other words,the set Λ(k) is such that: ##EQU1##

As shown above, the number of satellites attached to the node φ_(k) isdenoted n_(k).

The orbit concept discussed in the articles [6] and [7] has not yet beensubjected to any particular encoding process; an important feature ofthe invention is its proposed method of encoding Λ(k).

Neighborhood N(k)

Let φ_(k) be the k^(th) node of a surface S. The neighborhood N(k) ofφ_(k) is the subset N(k) of nodes equal to the union of φ_(k) and itssatellites. In other words:

    N(k)={φ.sub.k }∪Λ(k)

A neighbor of node φ_(k) is any node belonging to the set N(k). Thus:

    φ.sub.j εN(k) {φ.sub.j =Neighbor of φ.sub.k }

The neighborhood concept discussed in the articles [7] and [8] has notyet been subjected to any particular encoding process; an importantfeature of the invention is its proposed method of encoding N(k).

Atom A(R) and associated kernel

Let φ_(k) be the k^(th) node of a surface S. The atom of kernel φ_(k)A(R) is the pair (φ_(k),Λ*(k)) composed of:

the node φ_(k)

the set Λ*(K) of the addresses or the codes {A*(j1), . . . ,A*(jn_(k))}allowing the set of the corresponding atoms {A(j1), . . . ,A(jn_(k))}whose kernels are the satellites of φ_(k) to be retrieved from a memory.

In other words: ##EQU2##

The atom concept discussed in the article [6] has not yet been subjectedto any particular encoding process; an important feature of theinvention is its proposed method of encoding A(k).

Constraints and types of constraints

A constraint attached to an atom A(k) is any condition concerning thelocation of the node φ_(k) corresponding to the kernel of this atom. Forexample, in geology, the following constraints are of prime importanceto fitting a surface to precise or imprecise data observed or measuredby a manual or geophysical technique:

The coordinates (φ_(k) ^(x),φ_(k) ^(y),φ_(k) ^(z)) of the node φ_(k) areknown exactly. Such a constraint is called a control node typeconstraint and may be used in geology to encode the intersection of themodeled surface with a well.

The coordinates (φ_(k) ^(x),φ_(k) ^(y),φ_(k) ^(z)) of the node φ_(k) arenot known exactly but with a given certainty factor. Such a constraintis called a fuzzy control node type constraint and may be used to encodethe approximate location of a node of the modeled surface measuredapproximately by a given process. For example,

seismic data in geology, and

data corresponding to sections of organs observed with a microscope orusing ultrasound techniques or with a scanner in biology are of thistype.

The vector (φ.sub.λ -φ_(k)) joining the node φ_(k) to another nodeφ.sub.λ is approximately known with a given certainty factor. Such aconstraint is called a fuzzy vectorial link type constraint and is usedin geology, for example, to encode the throw vector of a fault, thenodes φ_(k) and φ.sub.λ being on opposite sides of the fault (see FIG.2).

The vector V orthogonal to the tangent plane to the surface S at nodeφ_(k) is approximately known with a given certainty factor. Such aconstraint is called a fuzzy vector normal type constraint and is usedin geology, for example, to encode the plane tangent to layer interfacesmeasured by a manual or geophysical technique (sounding).

Etc.

The certainty factors mentioned above are positive numbers whose valuesare assumed to be proportional to the a priori confidence that can beattached to the corresponding data. Accounting for these fuzzyconstraints constitutes one of the characteristics of the invention andis described in detail in section 2 (see sections 2.5 and 2.7.6),infra).

Notation

To simplify, no distinction will be drawn between the number of any nodeof a surface and the node itself. This implies that:

Λ(k)=set of satellites {φ_(j1), . . . ,φ_(jn).sbsb.k } of φ_(k) or setof corresponding numbers {j1, . . . ,jn_(k) },

N(k)=set of neighbors {φ_(k),φ_(j1), . . . ,φ_(jn).sbsb.k } of φ_(k) orset of corresponding numbers {k,j1, . . . ,jn_(k) }.

1.1.4 State of the art

Classical methods

A complete description of the "state of the art" concerning the modelingof complex surfaces is given in reference [7]. This article explainsprecisely why the classical methods used in automatic mapping and inComputer Aided Design (CAD) are not relevant to the modeling of complexsurfaces, as might otherwise be expected. The basics of the DSI method(also known as the DSI/DSA method) which overcomes these difficultiesare briefly introduced at the end of reference [7].

The DSI method

This method is based on the notion of "roughness" of a surface Scomposed of triangular and/or polygonal facets. This notion of"roughness" R(φ|k) at node φ_(k) of a surface S is defined as follows(See ref. [7], [8] and appendix): ##EQU3##

The positive or null coefficients {v.sup.α (k)} occurring in thisdefinition are weighting coefficients chosen by the user. Among theinfinity of possible choices, one of the most interesting, called"harmonic weighting", consists of choosing: ##EQU4##

The local roughnesses so defined are then combined in a global roughnessR(φ): ##EQU5##

The positive or null coefficients μ(k) occurring in this definition areused to modulate the contribution of the local roughness (R(φ|k) to theglobal roughness R(φ). If a uniform weighting is used, for example:

    μ(k)=1 k

Among other possible choices for the coefficients μ(k), the followingweightings can be used, where m is a given positive constant: ##EQU6##

The principle of the DSI method consists in computing the location ofeach node φ_(k) of the surface S to model in order to render it assmooth as possible while respecting the data and constraints governingthe shape of the surface. The paper [8] is devoted to a mathematicalpresentation of the DSI method and proposes a technique for minimizingthe global roughness while respecting the data and constraints.

Interactive version Of the DSI method

The method proposed in [8] for minimizing the global roughness is notwell suited to interactive use on an electronic computer (ormicroprocessor), which is why section 2, infra proposes a new method onwhich the invention is based. As mentioned previously, the object of theinvention is to propose an encoding and fitting process for the surfaceS based on this new version of the DSI method.

1.2 Examples of images of objects encoded and fitted with the process

1.2.1 Geological example

Oil prospecting is generally performed indirectly in the sense that oneis looking for geological layers which could have acted as a trap foroil; for example, deformations of salt layers called "salt domes" aregenerally excellent oil traps. In order to locate the oil precisely, itis necessary to know the exact shape of the "trap" layers andunfortunately, in many cases, until now there was no practical processusable to model it efficiently. The salt dome case is generally regardedas the most complex and that is why to prove the efficiency of thisencoding and fitting process FIG. 3 shows an image of such a surfaceobtained with it.

1.2.2 Biological example

Medical data acquisition processes (scanner or ultrasound) aremathematically identical to seismic techniques used in oil prospecting.For example, the ultrasound technique is identical to seismic reflectionand scanners are based on tomography as in seismic tomography. This isone reason why this encoding and fitting process may be used to modelthe skin of biological organs. FIG. 4 shows a model of the brain of a 7millimeters long human embryo obtained using this encoding and fittingmethod. Note that until now organs this small (the brain is 1 millimeterlong) could only be seen in the form of microscopic sections.

1.3 Encoding process for a surface S

1.3.1 Memory and memory address concepts

The term "memory" denotes any electronic computing device able to recorddata, codes or partial results and to retrieve them on request for usein subsequent operations. Moreover, the "address" of a memory is anycode identifying the physical location of the memory inside thecomputer.

1.3.2 Encoding an atom A(k)

The following process or one of its variants may be used to encode andstore each atom A(k) in a set of memories:

Encoding process

As suggested in FIG. 5, a contiguous memory area beginning at addressA*(k) stores consecutively the following data relating to the atom A(k)corresponding to the node φ_(k) :

1. Kernel=

either the coordinates (φ_(k) ^(x),φ_(k) ^(y),φ_(k) ^(z)) of the nodeφ_(k),

or the address or the code or any other data allowing the coordinates(φ_(k) ^(x),φ_(k) ^(y),φ_(k) ^(z)) of the node φ_(k) to be retrieved.

2. nb₋₋ sat=

either the number n_(k) of satellites linked to the node φ_(k),

or the address or the code or any other data allowing the number n_(k)of satellites linked to the node φ_(k) to be retrieved.

3. Sat=

Address of a memory area containing the addresses or the access codes{A*(j1), . . . ,A*(jn_(k))} allowing the n_(k) atoms {A(j1), . . .,A(jn_(k))} whose kernels {φ_(j1), . . . ,φ_(jn).sbsb.k } are thesatellites of the node φ_(k) to be retrieved.

4. Const=

Encoded constraints (see section 1.1.3) attached to the node φ_(k). Thisencoding will be explained in detail in section 1.1.3.

5. Info=

Memory area of variable size that can contain complementary dataconcerning the atom A(k). This data is specific to each particularapplication; for instance, in oil prospecting, if A(k) is an atom of thesurface separating two adjacent layers C1 and C2, then it may beinteresting to store in the Info field:

a list of physical properties (porosities, seismic velocities, etc) ofthe layer C1 at node φ_(k) corresponding to the atom A(k),

a list of physical properties (porosities, seismic velocities, etc) ofthe layer C2 at node φ_(k) corresponding to the atom A(k),

a list of geological properties (geological facies, presence of oil,etc) of the layer C1 at node φ_(k) corresponding to the atom A(k),

a list of geological properties (geological facies, presence of oil,etc) of the layer C2 at node φ_(k) corresponding to the atom A(k),

etc.

Variant 1

The memory area in which are stored the addresses or the access codes{A*(j1), . . . ,A*(jn_(k))} for retrieving the n_(k) atoms {A(j1), . . .,A(jn_(k))} whose kernels are the satellites {φ_(j1), . . .,φ_(jn).sbsb.k } of the node φ_(k) may be structured in two ways:

1. As shown in FIG. 6, the first way consists in using an array composedof n_(k) consecutive memories containing the addresses {A*(j1), . . .,A*(jn_(k))}.

2. The second way (FIG. 7) consists in using n_(k) arrays eachcontaining two consecutive memories such that for the α^(th) array:

the first memory contains the address A*(j.sub.α) of the atomcorresponding to the α^(th) satellite of the atom A(k),

the second memory contains the address of the (α+1)^(th) array. If α isequal to the number n_(k) of satellites (nb₋₋ sat=n_(k)), then thissecond memory is either unused or contains a code specifying that thelast satellite of A(k) has been reached.

Variant 2

It is possible not to store the nb₋₋ sat field but in this case, if thelist of addresses or access codes {A*(j1), . . . ,A*(jn_(k))} iS codedas a single array (see variant 1), it is necessary to add a memory areaA*(n_(k) +1) to store a code indicating that the list Λ*(k) is finished.

Variant 3

The Sat field has to be made big enough to contain the addresses oraccess codes {A*(j1), . . . ,A*(jn_(k))}.

Variant 4

Many variants are possible depending on the size of the Info field andon how it is partitioned. However, it is obvious that all these variantshave no influence on the use of the atom A(k) in connection with thefitting algorithm based on the DSI method described in section 2, infra.

Variant 5

Other fields can be added to the atom A(k) described above. Howeverthese other fields have no influence on the use of the atom A(k) inconnection with the fitting algorithm based on the DSI method describedin section, infra.

Variant 6

For any variant, what is really important for an efficientimplementation of the DSI method is that the address or the access codeA*(k) of a given atom A(k) should be the most direct path to the atomswhose kernels are the satellites of the kernel of A(k).

The main object of the method of encoding A(k) in a memory describedabove is to allow direct access to the satellites.

1.3.3 Encoding the Constraints

Each constraint relative to an atom A(k) will be stored in a specificmemory area whose description follows:

Process for encoding a constraint

As shown in FIG. 8, the encoding of a particular constraint attached toa given atom is performed in a memory area made up of two fields,Constraint and Next. These two fields are used as follows:

Constraint=

Memory area containing:

either the data related to the encoded constraint (see below),

or the address or the access code of a memory area containing the datarelated to the encoded constraint.

The data related to the encoded constraint must include a codeidentifying the type of the encoded constraint (see section 1.1.3).

Next=

Memory area containing:

either a code indicating that the current constraint is the lastconstraint attached to the atom A(k),

or the address of a memory area containing the next constraint (seebelow).

Process for storing all constraints attached to an atom

As described in section 1.3.2, encoding an atom A(k) requires a Constfield that will be used to store:

either a code specifying that the atom A(k) has no constraint,

or the address of a memory area containing the first constraint relatedto the atom A(k).

If the address of the n^(th) constraint related to the atom A(k) isstored in the Next field of the (n-1)^(th) constraint (n>1), then, asshown in FIG. 9, this implies that all the constraints related to A(k)are attached to A(k). This encoding process allows new constraints to beadded or old constraints to be removed without modifying theorganization of the other data.

Examples of data relating to a given constraint

Section 1.1.3 gives three examples of constraint types which are ofprime importance for the encoding and fitting of complex surfacesencountered in biology and geology. For each of these types ofconstraints data has to be attached to the Constraint field previouslydescribed and shown in FIG. 9:

For a fuzzy control node type constraint, in addition to the constrainttype, it is necessary to store at least the following data:

the approximate coordinates (φ_(k) ^(x),φ_(k) ^(x),φ_(k) ^(z)) of thenode φ_(k) corresponding to the kernel of the atom A(k) to which theconstraint is attached,

the certainty factor (positive number) proportional to the confidencethat can be placed in the approximate coordinates (φ_(k) ^(x),φ_(k)^(y),φ_(k) ^(z)).

For a fuzzy vectorial link type constraint, in addition to the type, itis necessary to store at least the following data:

the address of the node φ.sub.λ linked to the kernel φ_(k) of the atomA(k) to which the constraint is attached,

the three approximate coordinates of the vector (φ.sub.λ -φ_(k)) joiningthe node φ_(k) to the node φ.sub.λ,

the certainty factor (positive number) proportional to the confidencethat can be placed in the approximate coordinates of the vector (φ.sub.λ-φ_(k)).

For a fuzzy vector normal type constraint, in addition to the type, itis necessary to store at least the following data:

the three approximate components of the vector V orthogonal to thesurface S at the node φ_(k) corresponding to the kernel of the atom A(k)to which the constraint is attached,

the certainty factor (positive number) proportional to the confidencethat can be placed in the approximate coordinates of the vector V.

The three constraint types are described in detail in sections 2.5 and2.7.6, infra.

Variant

Until now only one certainty factor common to the x, y and z componentsof each constraint has been introduced but up to three differentcertainty factors can be assigned to these three components.

1.3.4 Encoding a surface S

In order to optimize the programming of the variant of the DSI methodpresented in section 2, infra, a surface composed of triangular and/orpolygonal facets is encoded as a set of N atoms whose kernels (nodes)are the vertices of the triangles and/or polygons. The data and codesfor each of these atoms {A(1), . . . ,A(N)} may be stored:

either in an array (see FIG. 10),

or in the Atom field of memory areas linked by their Next field (seeFIG. 11 ). In this case, the Next field associated with the atom numberk must contain the address or the access code of the next atom numbered(k+1); the Next field associated to the last atom A(N) must contain acode indicating that the last atom has been reached.

The second solution is better if it must be possible to add and/orremove atoms easily but requires more memory than the first solution.

1.3.5 Variants

To simplify the description, the memory areas used to encode the datahave been split up into fields and each of these fields has been given:

a name,

a memory location.

One of ordinary skill in the art would understand that the encodingprocess described in this patent does not depend on the names or theorder of these fields inside the memory areas used. Also, it is alwayspossible to add new fields for particular applications without modifyingthe storage process.

1.4 Fitting a surface S

1.4.1 Introduction

Let S be a surface encoded according to the process described in section1.3. There will be described hereafter a process based on this encodingto fit optimally the coordinates of some nodes {φ₁, . . . ,φN} in orderto:

make S as smooth as possible in relation to the global roughnesscriterion of the DSI method (see references [7], [8], the appendix andsection 1.1.4),

make S respect as much as possible the data and constraints, howeveraccurately they may be known, concerning the shape of this surface (seesections 1.1.1 and 1.1.3).

For that purpose, we use a new variant of the DSI method described inthe appendix and designed to use the encoding process described insection 1.3. Taking into account the measured or observed data, thismethod fits the three coordinates (φ.sub.α^(x),φ.sub.α^(y),φ.sub.α^(z))of each free node φ.sub.α successively and using an iterative method andthe following formula:

new value of φ.sub.α^(x) =ƒ.sub.α^(x) (old values of φ₁ ^(x), . . .,φ_(N) ^(x)),

new value of φ.sub.α^(y) =ƒ.sub.α^(y) (old values of φ₁ ^(y), . . .,φ_(N) ^(y)),

new value of φ.sub.α^(z) =ƒ.sub.α^(z) (old values of φ₁ ^(z), . . .,φ_(N) ^(z)).

Referring to the appendix describing the new version of the DSI method,these updating functions ƒ.sub.α^(x) (), ƒ.sub.α^(y) () and ƒ.sub.α^(z)() are derived from the local form of the DSI equations at point α andhave the following form (see sections 2.4.1, 2.4.2 and 2.7.3, infra):##STR1##

The nature and the role of the weighting coefficients {v.sup.α (k)} andμ(k) have been discussed in section 1.1.4 and are explained in article[8] and in section 2, infra. The values of the terms {Γ_(i) ^(x)α,Γ_(i)^(y)α,Γ_(i) ^(z)α }, {γ_(i) ^(x)α,γ_(i) ^(y)α,γ_(i) ^(z)α } and {Q_(i)^(x)α,Q_(i) ^(y)α,Q_(i) ^(z)α } depend on the types and the values ofthe constraints attached to the node φ.sub.α and their exact formulationis given in section 2, infra.

When the weighting coefficients {v.sup.α (k)} correspond to the harmonicweighting described in section 1.1.4 and in section 2, infra, theweighting functions ƒ.sub.α^(x) (), ƒ.sub.α^(y) () and ƒ.sub.α^(z) ()may be simplified and take the following form: ##STR2##

Hereafter there is described a method of calculating these functionsbased on the encoding described in section 1.3.

1.4.2 Notation

In order to simplify the description of the computation process thefollowing notation and definitions are used:

A block of operations is any list of operations beginning with an opencurly brace "{" and ending with a closing curly brace "}". Theoperations in a block may be individual operations or blocks ofoperations.

The term working memory or working variable refers to any memory areaused for storing intermediary results during a computation.

Loading the expression e in a memory m is the operation of evaluatingthe value of e and storing the result in the memory m. Such an operationis written as follows:

    m←e

For any memory m used to store a numerical value, incrementing m by theincrement i is the operation of computing the value m+i and storing theresult in the memory m. Such an operation is written as follows:

    m←m+i

In the above operation, the increment i may be:

a constant,

or the content of a memory,

or an arithmetical expression.

In the above operation, the increment i may be:

a constant,

or the content of a memory,

or an arithmetical expression.

1.4.3 Fitting process with harmonic weightings

When the coefficients {v.sup.α (k)} are the harmonic weightingsdescribed in section 1.1.4 and in the appendix, the values ƒ.sub.α^(x),ƒ.sub.α^(y) and ƒ.sub.α^(z) taken by the functions ƒ.sub.α^(x) (),ƒ.sub.α^(y) () and ƒ.sub.α^(z) () relating to the node φ.sub.α may becomputed efficiently by the following method based on the encoding of S.

0) Allocate the following working memories: ##EQU7##

1) Perform the following initializations (the order does not matter):##EQU8##

2) Let n.sub.α be the number of satellites of A(α). This number can beobtained from the nb₋₋ sat field of A(α).

3) Let Λ(α) be the set of the satellites of the atom A(α) accessiblefrom the Sat field of the atom A(α). For each satellite (kεΛ(α), repeatthe operations of the following block

{

3.1) Determine the number n_(k) of satellites of A(k). This number canbe obtained from the nb₋₋ sat field of A(k).

3.2) Determine the coordinates (φ_(k) ^(x),φ_(k) ^(y),φ_(k) ^(z)) of thenode φ_(k) of the atom A(k). These coordinates can be obtained from theNode field of A(k).

3.3) From the data in A(k) calculate or determine the value μ(k) of theweighting coefficient (see section 1.1.4 and section 2, infra).

3.4) Perform the following operations (in any order): ##EQU9##

3.5) Let Λ(k) be the set of satellites of the atom A(k) accessible fromthe Sat field of the atom A(k). For each satellite (βεΛ(k), β≠α) performthe operations of the following block:

{

3.5.2) Determine the current coordinates(φ.sub.β^(x),φ.sub.β^(y),φ.sub.β^(z)) of the kernel φ.sub.β of the atomA(β). These coordinates can be obtained from the Node field of A(β).

3.5.3) Perform the following three operations (in any order): ##EQU10##{

3.6) Perform the following three operations (in any order): ##EQU11## }

4) Perform the following operations: ##EQU12##

5) Let C(α) be the set of constraints attached to the atom A(α) that canbe accessed by the Const field of this atom. For each constraint(cεC(α)) repeat the operations of the following block:

{

5.1) Determine the type of the constraint c.

5.2) Depending on the type of the constraint c, perform the followingtwo operations:

5.2.1) Extract the data specific to the constraint c.

5.2.2) Depending on the data selected in (5.2.1), compute the values ofthe following expressions (see sections 2.5 and 2.7, infra): ##EQU13##

5.3) Perform the following operations (in any order): ##EQU14## }

6) Perform the following operations (in any order): ##EQU15##

At the end of this process, the working memories ƒ.sub.α^(x),ƒ.sub.α^(y) and ƒ.sub.α^(z) contain the values of the functionsƒ.sub.α^(x) (), ƒ.sub.α^(y) () and ƒ.sub.α^(z) () used to fit the valuesof the coordinates of the node φ.sub.α : ##EQU16##

This fitting process must be applied only to the "free" nodes, that isto say nodes that do not correspond to control nodes (see section1.1.3).

1.4.4 Variants

Variant 1

The implementation of the fitting method based on the DSI methoddescribed in the previous section in the case of harmonic weightingcoefficients may be very easily adapted to the case of any other kind ofweighting. The general structure of the process remains unchanged, theonly changes concerning the incrementing of the working memoriesallocated at step (0). How these working memories are incrementeddepends on the weighting adopted.

For example, if the following weighting coefficients {v.sup.α (k)} arechosen: ##EQU17## then it is easy to verify that the associated fittingfunctions are: ##STR3##

The process used to compute the values ƒ.sub.α^(x), ƒ.sub.α^(y) andƒ.sub.α^(z) of these functions at node φ.sub.α is then almost identicalto the process described for harmonic weighting:

0) Allocate the following working memories: ##EQU18##

1) Perform the following initializations (the order does not matter):##EQU19##

2) Determine the number n.sub.α of satellites of A(α). This number canbe obtained from the nb₋₋ sat field of A(α)

3) Let Λ(α) be the set of satellites of the atom A(α) accessible fromthe Sat field of the atom A(α). For each satellite (kεΛ(α)) repeat theoperations of the following block:

{

3.1) Determine the number n_(k) of satellites of A(k). This number canbe obtained from the nb₋₋ sat field of A(k). Then perform the followingoperation:

    n.sub.k.sup.2 ←n.sub.k ×n.sub.k

3.2) Determine the coordinates (φ_(k) ^(x),φ_(k) ^(y),φ_(k) ^(z)) of theNode φ_(k) of the atom A(k). These coordinates can be obtained from theNode field of A(k).

3.3) From the data in A(k) calculate or determine the value μ(k) of theweighting coefficient (see section 1.1.4 and section 2, infra).

3.4) Perform the following operations (in any order): ##EQU20##

3.5) Let Λ(k) be the set of satellites of the atom A(k) accessible fromthe Sat field of the atom A(k). For each satellite βεΛ(k), β≠α performthe operations of the following block:

{

3.5.2) Determine the current coordinates(φ.sub.β^(x),φ.sub.β^(y),φ.sub.β^(z)) of the kernel φ.sub.β of the atomA(β). These coordinates can be obtained from the Node field of A(β).

3.5.3) Perform the following three operations (in any order): ##EQU21##

}

3.6) Perform the following three operations (in any order): ##EQU22## }

4) Perform the following operations: ##EQU23##

5) Lot C(α) be the set of constraints attached to the atom A(α) accessedby the Const field of this atom. For each constraint cεC(α) repeat theoperations of the following block:

{

5.1) Determine the type of the constraint c.

5.2) Depending on the type of the constraint c, perform the followingtwo operations:

5.2.1) Extract the data specific to the constraint c.

5.2.2) Depending on the data selected in (5.2.1), compute the values ofthe following expressions (see sections 2.5 and 2.7, infra): ##EQU24##

5.3) Perform the following operations (in any order): ##EQU25## }

6) Perform the following operations (in any order): ##EQU26##

The numerical parameters stored in the Info field of the atoms may beinterpolated by the DSI method (see appendix: sections 2.1.2 and 2.2,infra). In this case, the process is almost identical to the processdescribed above used to fit a surface; the only modification consists inusing only one of the three functions normally dedicated to fitting thethree coordinates of the nodes and to calculate it according to theinterpolated parameter.

The present invention finds applications in the following fields inparticular:

in geology and in geophysics, it can be used to create a surface model,for example to construct a solid three-dimensional model, on which ageophysical simulation or an underground resources exploitationsimulation can be carried out; it can also be used to obtain a set ofdata for initializing a digital or analogue simulator of geophysical orresource exploitation phenomena; it can further be used, in conjuctionwith a graphics workstation, to display on a screen or to print out onpaper images used to monitor simulations as above;

in biology or in medicine, the invention can be used to constructthree-dimensional models of organs and/or prostheses or to simulate theeffects of plastic surgery; similarly, displayed or printedtwo-dimensional images can be used to monitor the above operations.

However, it is obvious that the invention applies more generally tomodeling surfaces of any kind and in particular any natural or man-madesurface. It can be of particular benefit in computer-aided design (CAD).

Finally, the present invention is in no way limited to the embodimentdescribed above and shown in the drawings and those skilled in the artcan vary and modify the embodiment as described above without departingfrom the scope of the invention.

NEW VERSION OF THE DSI METHOD FOR INTERACTIVE MODELING OF COMPLEXSURFACES

In Geology and Biology, a common problem is modeling complex surfacessuch as interfaces between areas of different kinds or having differentproperties. Classical modeling techniques based on Bezier interpolationsand spline functions [9] are not well suited to processing this type ofheterogeneous data. A different approach is based on the DSI ("DiscreteSmooth Interpolation") method [8]. In this approach, surfaces aremodeled using irregular triangular facets whose vertices must bedetermined to allow for the largest possible quantities of heterogeneousdata.

FIG. 12 shows an apparatus according to the present invention, which canbe utilized to implement the different surface modelling methods of thepresent invention. The surface modelling apparatus, as exemplified inFIG. 12, comprises a measuring device which measures to obtain a set ofgeometrical data, a meshing device which meshes a surface to bemodelled, a memory to store data related to and for each mode of themesh, a calculating device which processes fitting of the surfaceaccording to the stored data, and a surface representation generator forcreating a representation of the surface from fitted coordinates of eachnode.

2.1 Introduction: Statement of the problem

2.1.1 Surface S and related definitions

Referring to FIG. 1, let S be a surface composed of contiguous flattriangular facets embedded in the 3D Euclidian space (O|x,y,z). Thissurface need not be connex or closed. The following notation is used:##EQU27##

G is a graph whose nodes are the points of Ω and this set Ω isidentified with the set of integer indices {1, . . . ,N} used to numberthe nodes. The neighborhood N(k) of a node kεΩ iS defined as: ##EQU28##

In the following, it is assumed that these neighborhoods N(k) have thefollowing symmetry property:

    αεN(β) βεN(α)

The neighborhoods N(k) so defined thus generate a topology on G that canbe considered as an approximate topology for S itself.

2.1.2 Problem 1: Estimating a function φ defined on Ω

Assume that the position of each node kεΩ iS known and let φ(k) be afunction defined for all nodes k εΩ and known only on a subset of Ω:##EQU29##

Generally, L is different from Ω and there is an infinity of functionsφ(k) defined on Ω and interpolating the values {φ(I)=φ_(I) :IεL}. Thegoal is to select among all such functions the one which minimizes agiven criterion R*(φ) measuring:

the global roughness of the function φ(.),

the discrepancy between the function φ(.) and some constraints and datathat it should satisfy.

2.1.3 Problem 2: Fitting the surface S

In this case, only a subset of L of nodes kεΩ have a known positionφ(k)=φ_(k). ##EQU30##

As for the first problem, L is generally different from Ω and there isan infinity of surfaces S interpolating the points {φ(I)=φ_(I) :IεL}.The goal is to select among all such surfaces the one which minimizes agiven criterion R*(φ) measuring:

the global roughness of the surface S,

the discrepancy between the surface S and some constraints and data thatit must satisfy.

Allowing for the independence of the components (φ^(x) (k),φ^(y)(k),φ^(z) (k)) of each vector φ(k), it is easy to use the solution ofthe first problem to solve the second by simply writing:

    R*(φ)=R*(φ.sup.x)+R(φ.sup.y)+R*(φ.sup.z)

In the very particular case where the set Ω corresponds to the nodes ofa regular rectangular grid laid on the surface, consideration might begiven to using a Bezier, spline or Briggs type interpolation method (See[9], [1] and [3]); unfortunately these methods are not appropriate forirregular grids and cannot allow for the constraints described insections 2.2.3, 2.5 and 2.7.2.

2.2 Interpolating a function φ(k) defined on Ω

This section briefly introduces the Discrete Smooth Interpolation methoddescribed in [8].

2.2.1 Notation

In the following, φ denotes a column matrix of size N as: ##EQU31##

The problem does not depend on the method of numbering the nodes of thegrid. In order to simplify the notation, it will therefore be assumedthat a permutation of the elements of the matrix φ has been performed insuch a way that φ can be split into two submatrices φ_(I) and φ_(L) suchthat: ##EQU32##

Moreover, ∥ . ∥D denotes a seminorm associated to a square (N×N)positive semidefinite matrix [D] in such a way that, for any columnmatrix X of size N:

    ∥X∥.sub.D.sup.2 =X.sup.t ·[D]·X

2.2.2 Defining a global roughness criterion R(φ)

Let R(φ|k) be a local roughness criterion defined by the followingformula where {v.sup.α (k)} are given positive, negative or nullweighting coefficients: ##EQU33##

R(φ|k) can be used to derive a global roughness criterion R(φ) definedas follows, where μ(k) is a given weighting function defined on Ω:##EQU34##

R(φ) is completely defined by the coefficients {v.sup.α (k)} and μ(k)and it is easy to verify (See [8]) that there is always an (N×N)symmetric positive semidefinite matrix [W] such that:

    R(φ|k)=φ.sup.t ·[W(k)]·φ

2.2.3 Allowing for linear constraints in a least squares sense

Consider a square (N×N) matrix [A_(i) ] and an N column matrix B_(i)defining the linear constraint:

    [A.sub.i ]·φ≃B.sub.i

For an N×N positive semidefinite matrix [D_(i) ] "≃" means:

    ∥[A.sub.i ]·φ-B.sub.i ∥.sub.D.sbsb.i.sup.2 as small as possible

If there are several conditions of this type, the degree of violation ofthese constraints can be measured by the criterion ρ(φ) with: ##EQU35##2.2.4 Solution of the problem

Among all the functions φ(k) defined on Ω and interpolating the data{φ(I)=φ_(I) :IεL }, one is selected which minimizes the followingcriterion:

    R*(φ)=R(φ)+ρ(φ)

Developing R*(φ): ##EQU36##

The partition of the matrix φ induces a similar partition of thematrices [W*] and Q used to define R*(φ): ##EQU37##

The condition ∂R*(φ)/∂₁₀₀ .sbsb.I =[0] yields the following "DSIequation" characterizing all the functions {φ(k):kεΩ} constituting asolution of the problem: ##EQU38##

2.3 Uniqueness of the solution of the DSI equation

Definition

The set L of nodes where φ is known is said to be consistent relative toR*(φ) if each connex component of the graph G contains at least one nodebelonging to L.

Theorem

If the global roughness criterion R(φ) is such that: ##EQU39## the DSIequation based on R(φ) has a unique solution.

Proof

The theorem explained above is of theoretical interest only and does notaffect in any way the method described in this patent. It thereforerequires no proof here.

2.4 Local form of the DSI equation

In the following, instead of solving directly the DSI matrix equation,an iterative approach avoids the computation and storage of [W*]. Thisiterative approach is that used in this patent.

In order to simplify the notation, it will be assumed in the followingthat the positive semi-definite matrices [Di] of the DSI equation arediagonal. Moreover, when linear constraints of the type [A_(i) ]·≃B_(i)exist, the following notation is used: ##EQU40## 2.4.1 Computing∂R*(φ)/∂.sub.φ.sbsb.α

The definition of R(φ|k) implies: ##EQU41##

From this it can be deduced that: ##EQU42##

Allowing for the diagonal structure of [Di], it is easy to verify thatthe α^(th) element of ∂∥[A_(i) ]·φ-B_(i) ∥_(D).sbsb.i² /∂.sub.φ.sbsb.αis such that: ##EQU43##

Given the definition of R*(φ), it is possible to deduce that: ##EQU44##

The solution φ is such that ∂R*(φ)/∂.sub.φ.sbsb.α =0, hence the α^(th)component φ.sub.α of φ must satisfy the following equation: ##EQU45##

In this disclosure this equation is called the local form of the DSIequation at node α.

2.4.2 Proposition for an iterative algorithm

The above local form of the DSI equation suggests a straightforwardalgorithm to estimate the solution φ. For example, at iteration (n+1)the α^(th) component φ.sub.α.sup.(n+1) of the solution φ.sup.(n+1) mustsatisfy the DSI(α) equation so that an iterative algorithm can be:

let I be the set of indexes of nodes where φ.sub.α is unknown,

let φ be an initial approximate solution,

while (more iterations are needed). ##EQU46##

Note that this very simple algorithm does not use explicitly the matrix[W*] occurring in the DSI equation but has to compute repeatedlyproducts such as {v.sup.α (k)·v.sup.β (k)} which are the products usedto derive [W].

In fact, if the initial approximate solution is close to the exactsolution, few iterations are needed and the computation overhead becomesnegligible. For example, this occurs in an interactive context where theinitial solution φ is taken equal to the solution before some localmodifications are made by the user. In this case, despite the slightincrease in overhead, the local form of the DSI equation is preferableas it is much easier to use than the global form.

In order to show that the proposed algorithm actually converges, notethat R*(φ) can be expressed as a function of φ.sub.α :

    R*(φ)=A·φ.sub.α.sup.2 +B·φ.sub.α +C

The coefficients A, B and C are independent of φ.sub.α and: ##EQU47##

The minimum of the global criterion R*(φ) is achieved for the valueφ.sub.α =-B/(2A) which is precisely the value given by the DSI(α)equation. This shows that the algorithm converges because, at each stepof the iterative process, the value of the positive or null functionR*(φ) decreases.

2.5 Allowing for fuzzy data

This section introduces two types of linear constraints that will be ofspecial importance for Geometric Modeling, as will be shown in section2.7.6. The purpose of this section is to account for fuzzy dataconcerning the unknown values {φ.sub.λ :λεI}. Such data is assumed to beassociated with certainty factors ω² εR⁺ which are actually positiveweighting coefficients.

2.5.1 Case of isolated fuzzy data

Let λεI be a node index for which the unknown value φλ is assumed to beclose to an uncertain datum φλ with a certainty factor equal toω.sub.λ².

    φ.sub.λ ≃φ.sub.λ

where ≃ represents the following condition:φλ:

    ω.sub.λ.sup.2 ·|φ.sub.λ.sup.* -φ.sub.λ |.sup.2 minimum

Such a condition is easily accounted for in the DSI method as the i^(th)constraint with: ##EQU48##

The corresponding coefficients Γ_(i).sup.α, Q_(i).sup.α and γ_(i).sup.αused in the local DSI(α) equation are thus defined by: ##EQU49## 2.5.2Case of differential fuzzy data

Let λεI be a node index for which the unknown value φλ is assumed to belinked to another value φμ corresponding to another index μ≠λ. Assumethat this link is of the following form, where Δ₈₀ μ is a given value:

    (φ.sub.μ -φ.sub.λ)≃Δ.sub.λμ

This relation is assumed to be fuzzy and is modeled through thefollowing condition involving a certainty factor ω.sub.λμ² :

    ω.sub.λμ.sup.2 ·|(φ.sub.μ.sup.* -φ.sub.λ.sup.*)-Δ.sub.λμ |.sup.2 minimum

Such a condition can easily be accounted for in the DSI method as thei^(th) constraint with: ##EQU50##

The corresponding coefficients Γ_(i).sup.α, Q_(i).sup.α and γ_(i).sup.αused in the DSI(α) equation are thus defined by: ##EQU51## 2.5.3Choosing the certainty factors

Consider the term M.sub.α^(*) of the DSI(α) equation: ##EQU52##

In the two examples discussed above, the terms γ₁.sup.α are either equalto zero or equal to the certainty factor:

    γ.sub.i.sup.α =ω.sub.i.sup.2

This suggests choosing for ω_(i) ² a given percentage P_(i) of M.sub.α :

    ω.sub.i.sup.2 =p.sub.i ·M.sub.α p.sub.i >0

The term M.sub.α^(*) of the DSI(α) equation is then:

    M.sub.α.sup.* =M.sub.α ·(1+p.sub.1 (α)+p.sub.2 (α)+ . . . )

where p_(i) (α) is either equal to a given percentage P_(i) or equal tozero.

2.6 Choosing the weighting coefficients

The choice of the weighting coefficients {v.sup.α (k)} and {μ(k)} iscompletely free except that the {μ(k)} coefficients have to be positiveor equal to zero. In the following an example is given of how to choosethese two families of coefficients.

2.6.1 Choosing the {v.sup.α (k)} harmonic weighting coefficients

Let Λ(k) be the "orbit" of k defined as follows:

    Λ(k)=N(k)-{k}

Let |Λ(k)| be the number of elements of (k). The weighting is calledharmonic weighting if the coefficients {v.sup.α (k)} are chosenaccording to the following definition: ##EQU53##

Harmonic functions have the characteristic property that they are equalat any point to their own mean on a circle centered on this point. Asone can see, R(φ|k)=0 if φ_(k) iS equal to the mean of the valuesφ.sub.α surrounding the node k, this is why, by analogy with harmonicfunctions, the term harmonic has been adopted for the weightingcoefficients {v.sup.α (k)} described in this section.

If the DSI(α) equation is developed with these coefficients, thefollowing equation is obtained, which can be easily translated intoprogramming language: ##EQU54## 2.6.2 Choosing the {μ(k)} coefficients

The coefficients {μ(k)} have been introduced in order to modulatelocally the smoothness of the interpolation. If there is no specialreason to do differently, a uniform weighting may be used:

    μ(k)=1 kεΩ

Among all non-uniform weighting schemes, one seems to be particularlyinteresting if smoothness is needed in the neighborhood of the data:##EQU55##

2.7 Application to Geometric Modeling

This section shows how to use the DSI method to estimate thetriangulated surface S itself for which it will be assumed that partialdata is available composed of:

the exact location of some apices of the triangles,

the approximate location of some apices of the triangles,

vectorial constraints specifying the shape of S,

etc.

To this end, φ.sub.α is defined as the current vector joining the originof R³ to the current node of Ω and φ^(x), φ^(y) and φ^(z) denote itscomponents on the {ox,oy,oz} orthonormal basis of R³ :

    φ=(φ.sup.x,φ.sup.y,φ.sup.z)

2.7.1 Defining a global roughness criterion R(φ)

As usual, the set of nodes Ω is split into two subsets I and L suchthat: ##EQU56##

The nodes IεL whose location φ_(I) εR³ is known will be called controlnodes and the aim is to determine the location of the remaining pointsiεI in such a way that the following local criterion R(φ|k) based on agiven set of weighting coefficients {v.sup.α (k)} is as small aspossible: ##EQU57##

This criterion constitutes the vectorial form of the local DSIcriterion. The associated global vectorial DSI criterion R(φ) is definedas follows where {μ(k)} are given non-negative weighting coefficients:

    R(φ)=Σμ(k)·R(φ|k)

Defining R(φ^(x)), R(φ^(y)) and R(φ^(z)) as previously: ##EQU58##

Pythagoras' theorem yields:

    R(φ)=R(φ.sup.x)+R(φ.sup.y)+R(φ.sup.z)

2.7.2 Allowing for linear constraints in a least squares sense

Consider the following matrices sized so that N is the total number ofnodes of the set Ω: ##EQU59##

To have the components φ^(x),φ^(y) and φ^(z) satisfy the followingconstraints

    [A.sub.i.sbsb.x.sup.x ]·φ.sup.x ≃B.sub.i.sbsb.x.sup.x

    [A.sub.i.sbsb.y.sup.y ]·φ.sup.y ≃B.sub.i.sbsb.y.sup.y

    [A.sub.i.sbsb.z.sup.z ]·φ.sup.z ≃B.sub.i.sbsb.z.sup.z

the following criterion is introduced: ##EQU60## 2.7.3 Solution of theproblem

Among all the surfaces S interpolating the data {φ(I)=φ_(I) :IεL }, theaim is to select the one minimizing the following criterion:

    R*(φ)=R(φ)+ρ(φ)

It is easy to verify that this equation can be rewritten as follows:##EQU61##

Allowing for independence of the components φ^(x),φ^(y) and φ^(z) :##EQU62##

Minimizing the vectorial DSI criterion is thus equivalent to minimizingthe corresponding three DSI criteria applied to the three components ofthe current vector φ.

2.7.4 Note

If the surface S is composed of several disjoint patches the theoremexplained in section 2.3 ensures that there exists a unique solution tothe problem if at least one triangle apex within each patch has a knownposition.

2.7.5 Interactive use

In Geometric Modeling, the goal is to define interactively the geometricshape of the graph G associated with the set of nodes Ω. To this end, itis assumed that an initial shape is known and the user moves some nodesand changes interactively some constraints (control nodes and fuzzydata); these modifications are then propagated to the nodes iεI throughthe vectorial form of the DSI method.

As explained previously, this problem is broken down into threesub-problems corresponding to the three coordinates (x,y,z) and the bestway to solve these three problems is to use the iterative algorithmassociated with the local form of the corresponding DSI equations.

Between steps of the modeling process, not only can the position of thecontrol nodes be modified but also the subset L itself can be changed ifnecessary; for example, to apply DSI only to a subset Γ of nodes, it isnecessary to choose L in such a way that Γ Γ where Γ stands for thecomplementary set of Γ in Ω.

Moreover, because the problem is broken down into sub-problems solvedindependently and corresponding to the three coordinates, differentsubsets L_(x), L_(y) and L_(z) can be defined for each coordinate ifnecessary; for example, if the (x,y) coordinates must not be modified itis sufficient for L_(x) and L_(y) to be equal to the whole set Ω.

2.7.6 Allowing for fuzzy geometrical data

Initially, the DSI method was designed to model complex surfacesencountered in the field of natural sciences, for example in geology andbiology. In this case, there is generally much imprecise data to takeinto account such as: ##EQU63##

Fuzzy control nodes

By definition, a "fuzzy control node" is any node λεI whose position φλis known as being approximately equal to a given vector φλ with a givendegree of certainty ω.sub.λ². Referring to section 2.5.1, such fuzzydata can be taken into account with the DSI method by introducing thefollowing constraints in the R*(φ^(x)), R*(φ^(y)) and R*(φ^(z))criteria: ##EQU64##

Fuzzy vectorial constraints

In many applications it is necessary to control the vector (φ.sub.μ-φ.sub.λ) joining two nodes φ.sub.λ and φ.sub.μ of a surface. Bydefinition, such a constraint is called a "fuzzy vectorial constraint"if it must be satisfied with a given degree of certainty ω.sub.λμ². Someexamples of constraints of this type are shown in FIG. 2.

In practice, at least one of the two points φ.sub.λ and φ.sub.μ has anunknown position, say the one corresponding to λεI, and two importantcases have to be considered:

The first case corresponds to FIGS. 2a and 2b where the shape of anobject must be controlled in such a way that (φ.sub.μ -φ.sub.λ) is equalto a given vector Δ.sub.λμ.

The second case corresponds to FIG. 2c where the shape of an object mustbe controlled in such a way that (φ.sub.μ -φ.sub.λ) is orthogonal to agiven vector V. This situation occurs when the vector normal to asurface is to be controlled; in this case, if V is the normal vector atnode λ of the surface, it is wise to control its orthogonality with allthe vectors (φ.sub.μ -φ.sub.λ) for all μεN(λ).

Referring to section 2.5.2, the first case is easily taken into accountin the DSI method by introducing the following constraints in theR*(φ^(x)), R*(φ^(y)) and R*(φ^(z)) criteria with a weighting coefficientequal to ω.sub.λμ² : ##EQU65##

For the second case, the problem is not very different because:##EQU66##

In this case, the values (φ.sub.λ^(x),φ.sub.λ^(y),φ.sub.λ.sup.z) and(φ.sub.μ^(x),φ.sub.μ^(y),φ.sub.μ^(z)) on the righthand side of theequations can be set to the corresponding values inferred at theprevious step of the iterative process; if

    Δ.sub.λμ.sup.x =-{V.sup.y ·(φ.sub.μ.sup.y -φ.sub.λ.sup.y)+V.sup.z ·(φ.sub.μ.sup.z -φ.sub.λ.sup.z)}/V.sup.x

    Δ.sub.λμ.sup.y =-{V.sup.z ·(φ.sub.μ.sup.z -φ.sub.λ.sup.z)+V.sup.x ·(φ.sub.μ.sup.x -φ.sub.λ.sup.x)}/V.sup.y

    Δ.sub.λμ.sup.z =-{V.sup.x ·(φ.sub.μ.sup.x -φ.sub.λ.sup.x)+V.sup.y ·(φ.sub.μ.sup.y -φ.sub.λ.sup.y)}/V.sup.z

the second case becomes identical to the first case.

Note

In allowing for fuzzy data as described above only one certainty factorω.sub.λμ² is used for all three components of the constraints. Inpractice this is not mandatory and, because of the independence of thethree components, it is always possible to use three different certaintyfactors ω.sub.λμ^(x2), ω.sub.λμ^(y2) and ω.sub.λμ^(z2) corresponding tothe three components.

2.8 Generalization

Before concluding, some straightforward generalizations of the approachare worth mentioning:

The DSI method can be used on any kind of polygonal facets. In practice,triangles are recommended since they are the simplest polygons and areeasy to handle in a computer program.

The DSI method can be used to estimate polygonal curves in a 3D space;in this case, triangles are replaced by segments.

2.9 Conclusion

The local form of the DSI equation provides a simple and powerful toolfor modeling complex surfaces encountered in Geology and Biology, forexample. In particular:

There are no restrictions as to the mesh used to model a surface and itis possible to uses automatic algorithms to build them; for example, ifthe mesh is composed of triangles, it is possible to use algorithmsderived from the Delaunay method (See [2] and [5]). Moreover, the sizeof the mesh can easily be fitted locally to the complexity of thesurface.

The DSI method can cope with large quantities of precise or fuzzy data.This may be particularly useful in Natural Sciences where the aim is notto generate aesthetic surfaces but to fit precise and fuzzy dataconsistently.

The algorithm derived from local DSI equation is fast enough to allowinteractive use.

Nothing has been said about the representation of the facets. Dependingon the application, plane facets may suffice. However, facets can ifnecessary be represented by non-planar surface patches interpolating thenodes of the mesh (See [4]).

The invention being thus described, it will be obvious that the same maybe varied in many ways. Such variations are not to be regarded as adeparture from the spirit and scope of the invention, and all suchmodifications as would be obvious to one skilled in the art are intendedto be included within the scope of the following claims.

I claim:
 1. A process for modeling a surface representing an interfacebetween two areas with different properties in a three-dimensional body,comprising the steps of:measuring to obtain a set of geometrical data,from said interface of said three-dimensional body, relating to thesurface and associated with respective points on said surface; meshingthe surface so that all said points are a subset of nodes of a mesh;storing at a specific memory address for each node of the mesh, thefollowing data:first, second, and third coordinates, a number ofsatellite nodes, satellite node address data providing access tospecific coordinates of said satellite nodes and consequently tosatellite node related data relating thereto, and geometrical dataassociated with said each node of the mesh; for each node of the mesh,defining a local roughness index derived from a weighted sum of thecoordinates of the node and of its satellites; defining a globalroughness index obtained by summing the local roughness indicesassociated with the nodes, and a global violation index of saidgeometrical data, to thereby generate a sum of said global roughnessindex and said global violation index; fitting fitted coordinates ofeach node processed by an iterative method in which for each step of aniteration there is added a weighted combination of the coordinates ofthe node and of the satellites of said node and a combination of thegeometrical data associated with said node, to thereby minimize said sumof the global roughness index and the global violation index; andcreating a representation of the surface from the fitted coordinates ofeach node.
 2. A process according to claim 1, wherein:the surface ismeshed using triangles.
 3. A process according to claim 1 wherein:thefitting step comprises three separate sub-steps respectively using thefirst, second, and third coordinates of the nodes and the coordinates ofthe satellite nodes corresponding thereto, respectively.
 4. A processaccording to claim 1, wherein:the geometrical data includes vector databetween two given nodes.
 5. A process according to claim 1, wherein:thegeometrical data includes data of a vector normal to the surface to bemodeled.
 6. A process according to any one of claims 4 to 5, wherein:atleast one geometrical datum is associated with a coefficientrepresenting a certainty factor by which said at least one datum isknown, the global violation index is a function of said certaintyfactors.
 7. A process according to claim 1, further comprising the stepsof measuring properties of the three-dimensional body in a region of atleast one node of the mesh and storing the measured properties ataddresses specific to the at least one node so as to obtaincomplementary data of the three-dimensional body.
 8. A process accordingto claim 1, wherein:for each local roughness index, a ratio between afirst and a second weight of a corresponding node equals a negative ofthe number of satellite nodes.
 9. A process according to claim 8,wherein:said weighted combination uses substantially the same weightingas the ratio.
 10. A process according to claim 1, wherein:the summing ofthe local roughness indices is a weighted summing.
 11. A processaccording to claim 10, wherein:said weighted combination uses theweighting associated with the local roughness indices.
 12. A processaccording to claim 1, wherein:the three-dimensional body is a geologicalformation.
 13. A process according to claim 1, wherein:the surfacerepresents variations of a property of a two-dimensional body in avicinity of a plane intersecting said body, said property beingrepresented by a remaining dimension, and a representation of thesurface is used to optimize visualization of underground resources ofsaid three-dimensional body, said three dimensional body being ageological formation.
 14. A process according to claim 1, wherein:therepresentation of the surface is a graphical representation on a flatmedium.
 15. A method according to claim 1, wherein:the representation ofthe surface is a three-dimensional model.
 16. A device for modeling asurface representing an interface between two areas with differentproperties in a three-dimensional body, comprising:means for measuringto obtain a set of geometrical data, from said interface of saidthree-dimensional body, relating to the surface and associated withrespective points on said surface; means for meshing the surface so thatall said points are a subset of nodes of a mesh; means for storing at aspecific memory address for each node of the mesh, the followingdata:first, second, and third coordinates, a number of satellite nodes,satellite node address data providing access to specific coordinates ofsaid satellite nodes and consequently to satellite node related datarelating thereto, and geometrical data associated with said each node ofthe mesh; calculation means for fitting fitted coordinates of each nodeprocessed by an iterative method in which for each step of an iterationthere is added a weighted combination of the coordinates of the node andof the satellites of said node and a combination of the geometrical dataassociated with said node, to thereby minimize a sum of a globalroughness index obtained by summing local roughness indices associatedwith the nodes, each local roughness index derived from a weighted sumof the coordinates of the node and the coordinates of its satellites,and of a global violation index of said geometrical data; and means forcreating a representation of the surface from the fitted coordinates ofeach node.